Attaching packages

knitr::opts_chunk$set(echo = TRUE, 
                      warning = FALSE, 
                      message = FALSE)
library(tidyverse) # Default.
library(here) # To create file paths to data within different folders.
library(sf) # To read in spatial data and make simple features.
library(raster) # To work with raster data.
library(fasterize) # To make simple features and geometry.

Reading in the raster data

NOTE: Working with raster data and images on the R Bren Server 3.6 will produce chunky, pixelated raster images on the .Rmd, but will produce clean, smooth raster images when knitted to an html.

knitr::include_graphics(here("img", 
                             "landsat.png"))

Loading in and inspecting the data

landsat_file <- here("data", 
                     "Landsat7.tif")

ls_1 <- raster(landsat_file)
ls_1
## class      : RasterLayer 
## band       : 1  (of  5  bands)
## dimensions : 1758, 3701, 6506358  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : -59564.57, 51465.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : /home/csimms/esm_244_adv_data_analysis/labs/lab_6/data/Landsat7.tif 
## names      : Landsat7 
## values     : 0, 255  (min, max)
plot(ls_1)

ls_2 <- raster(landsat_file, 
               band = 2)
ls_3 <- raster(landsat_file, 
               band = 3)
ls_4 <- raster(landsat_file, 
               band = 4)

ls_stack <- raster::stack(landsat_file) # To read all the layers in at once.
ls_stack
## class      : RasterStack 
## dimensions : 1758, 3701, 6506358, 5  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : -59564.57, 51465.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## names      : Landsat7.1, Landsat7.2, Landsat7.3, Landsat7.4, Landsat7.5 
## min values :          0,          0,          0,          0,          0 
## max values :        255,        255,        255,        255,        255

Preparing the data by overwrighting the original layers

ls_1 <- raster::aggregate(ls_1, 
                          fact = 3, 
                          fun = mean)
ls_2 <- raster::aggregate(ls_2, 
                          fact = 3, 
                          fun = mean)
ls_3 <- raster::aggregate(ls_3, 
                          fact = 3, 
                          fun = mean)
ls_4 <- raster::aggregate(ls_4, 
                          fact = 3, 
                          fun = mean)

plot(ls_1, 
     col = hcl.colors(n = 100, 
                      palette = 'Blues 2')) # To base a color scheme on a single palette.

plot(ls_2, 
     col = hcl.colors(n = 100, 
                      palette = 'Greens 2'))

plot(ls_3, 
     col = hcl.colors(n = 100, 
                      palette = 'Reds 2'))

plot(ls_4, 
     col = hcl.colors(n = 100, 
                      palette = 'Reds 2'))

sbc_rast <- raster(here("data", 
                        "county.tif"))

plot(sbc_rast)

plot(ls_3)

# ERROR Code in Lab Playthrough:

# mask(ls_3, sbc_rast) %>% plot()

# ls_3 <- mask(ls_3, sbc_rast)

# ls_4 <- mask(ls_4, sbc_rast)

Working with raster algebra

vec1 <- 1:5
vec1
## [1] 1 2 3 4 5
vec1*2
## [1]  2  4  6  8 10
vec1^2
## [1]  1  4  9 16 25
ls_3
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : Landsat7 
## values     : 8.444444, 255  (min, max)
ls_3*2
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : Landsat7 
## values     : 16.88889, 510  (min, max)
log(ls_3)
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 2.133509, 5.541264  (min, max)
plot(ls_3); plot(log(ls_3)) # To compare.

vec2 <- 6:10
vec1+vec2
## [1]  7  9 11 13 15
ls_3+ls_4
## class      : RasterLayer 
## dimensions : 586, 1234, 723124  (nrow, ncol, ncell)
## resolution : 90, 90  (x, y)
## extent     : -59564.57, 51495.43, -404675.9, -351935.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 35.33333, 498.7778  (min, max)

Using ‘raster::calc()’

ls_stack <- stack(ls_1, 
                  ls_2, 
                  ls_3, 
                  ls_4)

ls_mean <- raster::calc(ls_stack, 
                        fun = mean, 
                        na.rm = FALSE)

plot(ls_mean)

Analysis: NDVI

knitr::include_graphics(here("img", 
                             "spectrum.png"))

knitr::include_graphics(here("img", 
                             "ir_photo.jpg"))

NDVI equation

\[NDVI = \frac{NIR - Red}{NIR + Red}\]

ndvi <- (ls_4 - ls_3) / (ls_4 + ls_3)

plot(ndvi, 
     col = hcl.colors(100, 
                      palette = 'Grays'))

is_forest <- function(x, 
                      thresh = 0.3) {
  y <- ifelse(x >= thresh, 
              1, 
              NA)
  return(y)
}

forest <- calc(ndvi, 
               fun = is_forest)

plot(forest, 
     col = 'green')

‘ggplot’ and rasters

ndvi_df <- raster::rasterToPoints(ndvi) %>% 
  as.data.frame()

forest_df <- raster::rasterToPoints(forest) %>% 
  as.data.frame()

ggplot(data = ndvi_df, 
       aes(x = x, 
           y = y, 
           fill = layer)) + 
  geom_raster() + 
  geom_raster(data = forest_df, 
              fill = "green") + 
  coord_sf(expand = 0) + 
  scale_fill_gradient(low = "black", 
                      high = "white") + 
  theme_void() + 
  theme(panel.background = element_rect(fill = "slateblue4"))